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1. INTRODUCTION 


O A E, A AAA A nn 


Many familiar physical and operational processes are 
well described, in whole or part, by examining their "event 
streams" cver time. Some such well-known processes are: the 
plow, Of traffic through ап intersection; the arrival of 
persons at scme service facility such as a bank teller's 
window, a service station fuel pump or a grocery store check 
out counter; and, the arrival of telephone calls or radio 
transmissions at sone switchboard or other type of 
communications terminal.  Analogous processes abound bcth in 
nature and in the course of our everyday lives. Ey the 
proper definition of an "event" in these situations, the 
process being observed will be characterized by the 
probabilistic nature of the flow, over time, of the events 
of which it is composed. In the above three examples the 
events could be defined respectively as the arrival of a 
vehicle at the intersection, the arrival of a customer at 
the service facility, and the arrival of the telephone call 
Or radio transmission at the terminal. The process may then 
be analyzed ry examining the interaction of various event 
streams with different intersection configurations, service 


policies and terminal capacities. 


A commen method used to perform such analyses is to 
consider the event streams to be homogeneous. This could 
mean that the expected number of events to cccur in any two 
Or more time intervais of equal length is the same (simple 
Memo genelty), or that the distribution of the nunter of 
events occurring in any two or more equal time intervals is 
the same (complete homogeneity). The homogeneous Poisson 


process is often used as a tool in the analysis of such 


Ta 





activities. For very simple systems involving event streams 
the use of the homogeneous Poisson process as a model leads 
to tractable analytical results. Most systems of interest, 
however, are not amenable to purely analytical methods, and 
simulaticn of the processes by digital means becomes 


necessary. 


The assumption of homogeneity in event streams is often 
a very restrictive one. Тие rsh "hour" phenomenon, 
well-known and abundantly cursed by motorists, provides 
cogent evidence that the modeling of event streams is not 
always weil served by the homogeneity assumption. The 
intensity of event streams varies over time for many 
physical cr operational processes. In these cases, purely 
analytical methods must be abandoned almost immediately in 
favcr of simulation techniques. (The intensity of the event 
Stream is defined to be the derivative with respect to t of 
the expected number of events in an interval of length 
(0,t ]-) 


The ncn-hcmogeneous Poisson process is often employed in 
the analysis of processes that exhibit gross departures from 
the homogeneous event stream criterion. If these processes 
are to be simulated, it is necessary to first describe the 
nature of the inhomogeneity (i.e. how does the intensity of 
the event stream vary over time?) and then to artificially 
generate event streams that behave in accordance with the 


descripticn. 


Of course there are infinite variations in the tyres of 
inhomogeneities that can occur or that сап be construed. 
However it is intuitively appealing to consider event 
Srzeams that display varying intensities of the following 


types: 


i) increasing continuously over time; 
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ii) decreasing continuously over tine; 


iii) increasing and then decreasing continuously over 


time; 


iv) decreasing then increasing continuously over time. 


A ncn-hcmogeneous Poisson process that has quadratic 
properties can be manipulated to produce the above-mentioned 
effects. The effective simulation of such a process on the 
computer is the subject of this thesis and has been 
nctivated by the work of P. A. W. Lewis, Professor, Naval 
postgraduate School, апа 6. 5. Shedler, IBM Research 
шарогатогу. 


There are of course other types of inhomogeneities, such 
as cyclic variations cane Of day effect), but-we do not 
consider them here. They have been discussed by СОХ 
moc | папа LEWES ff Ref. 9]. Em particular LEWIS [Ref. 9] 
describes a process consisting of arrivals at an intensive 
се Unit in a hospital. It is shown empirically that, in 
Egon to the time of day cycle, long term fluctuaticns in 
the intensity function can be adequately described by an 


intensity function whose logarithm is quadratic. 


ШЕТ апа SHEDLER [ Ref. 11] proposed a new methed for 
generating a non-homogeneous Poisson process with an event 
Stream intensity (rate) function that is of degree-two 
exponential polynomial form. (The use of exponential 
polynomials is natural in this context since an intensity 
function is a positive function.) The new method aprears to 
have the virtue of increased efficiency over the more 
conventicnal time-scale transformation technique when 


implemented on a high speed digital computer. 


ВЕСЕ Су in chis "context is measured in terms of 
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computer memory requirements and computational speed. The 
prcblem of efficiency comparison is recognized by LEWIS and 
SHEDLER in the final pages of their paper (Ref. 11, p. 15]: 
"There remains the question of efficiency (cf the proposed 
algorithm) . . . for generation of a non-homogeneous Pcisson 
process with degree-two exponential polynomial rate 
function, relative to generation via tine scale 


transformation of a homogeneous Poisson process." 


After scme brief discussion of the requirements of 
igplementing the time-scale transformation algorithm the 
report concludes . . +. "We therefore would expect the exact 
method of (the provosed algorithm) to be much faster, 


although at the expense of some complexity of programming." 


The objective of this author has been: to implement 
ШЕП che algorithm of LEWIS and SHEDLER [ Ref. 11] and the 
conventicnal time-scale transformation algorithm on the IBM 
360/67 Ccmputer System in FORTRAN IV language; to define 
reasonable measures of effectiveness for comparing the two 
algorithms; and to determine which algorithm is the more 


КЫ леп їп terms OF the measures of effectiveness defined. 


Section II discusses the definition of a non-homogeneous 
Pcisson process. It also states some special properties of 
Pcisscn processes that are used in the development cf the 
algorithms investigated in this thesis, and concludes with a 
general discussion of two basic methods of generating a 


non-homogenecus Poisson process. Section III gives a step 


by step description of the two algorithms that were 
implemented into computer programs. The method ОЁ 
implementation is the topic of Section IV.  Methcds of 
comparing the algorithms are presented in Section M 


Section VI presents conclusions drawn from the comparison 
mma Makes a  reconmendation for improving the LEWIS and 


ESSDLER [Ref. 11) algorithm. ther recommendations for 
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further study are also listed in Section VI. Appendixes 
provide additional details that would have been awkward to 
include in the main body of the thesis. Computer program 


listings are provided after Appendix E. 
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II. DEFINITION AND PROPERTIES OF NON-HOMOGENEOUS PCISSON 


PROCESSES 
A. GENERAL 
ла = section will present basic definitions and 
explanations concerning the concept of non-homogeneous 
Pcisson processes. References cited may be consulted if a 
more in-depth understanding is desired. Only the 
fundamental concepts necessary for understanding the 


specific non-homogeneous Poisson process under consideration 


are presented. 


DEF INITICNA OF A NON-HOMOGENEOUS POISSON PROCESS 


LEWIS and SHEDLER [Ref. 11] define the non-homogeneous 


Poisson process on the real line as follows: 


1. The number of events in any finite set o 


non-overlapping intervals are independent random variables. 


De Let A (t) be a monotone non- decrsasing 
megnt-continucus function which is bounded or any finite 


interval. Then the number of events in any interval, e.g. 


(0,50), ћаз а Poisson distribution with paramerer 
A (tg) zm ae 
The function A(t) - A(0) is called, among other things, the 


l6 





mean value function since the expected number of events in 
Acera, E) 15 just A(t) - A(0). Property 1 is the 
independent ircrement property of the Poisson process; it is 
basic to the idea of a Poisson process. Figure 1 provides a 


graphical representation of the definition above. 


The above Черт оп insures that 
ЖОО) = (0)}/{ A (t9) = A(0)) meets the following 
Ens cerra for ап arbitrary function F(t), to be a valid 
Hcctra2butron function on the interval (0,54), Cran. LARSON 
NICE. 7, ch. 3]; 


MO EE) EL 


Дет) ШШЕ SO Lim PE TÉ) = 1, 0 < € < t 


m 0 
ELO t Es 


Ea ЕО) for all a< b іп (0,50) 
іу) IN (БОП) -СЕ(Р) for all р  (0,t 


pu 


where h > Q 


07” 


же - па 2(5) = СГ) (Р) - A(0)7^7 f A (to) О E 5011045 
muat if Is absolutely Continuous in (0,2) then 
EEUU /dt - A(t)/f A ($9) - A(0)) is a valid density function 
on tne interval (0, tQ]. Mic пе олљ a A(t is called che 


intensity function (or rat= function) of the process and 


hate) = x E{number of events in (0,5,1! 
md 
a ae {A(t) = A(0)]} 
d 
EN ut) 


п 





A(t) 
Alto) 
Alta) 


Alta) 


A(t) 


A(t, ) 
MO) ZT --- ----- -- > 2 22020 2 rn 2 


Situation: Events occur randomly in time (i.e. along t-axis) 
1) If X = number of events in fixed interval (t,,t,] 

Y = number of events in fixed interval (t,,t2] 

Z = number of events in fixed interval (t,,t,] 


S = number of events in fixed interval (0,4,1; 


Then X, Y and Z must be independent random variables (note: 


3 and 2 are not independent because of overlap in the intervals). 


2) Given the monotonic increasing right continuous function 


A(t), the number of events in any fixed interval (t; tj] 


must have the Poisson distribution with parameter A(t,) - A(t,). 


Th 
de X ^ Poisson(A(t,)) Е A(t, )} 
Y ^ Poisson(A(t.) - A(t,)) 
Z ~ Poisson{A(t,) - ACC Q)) 
S > Poisson{A(t,) - A(0)} 
Fiaure 1 - DEFINITION Or A NON-HOJOGENEOUS POISSON EROCESS 


(SOPA AD ZUPRESENTATION) 
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EEN TE NSTITY FUNCTION 


The mest intuitively appealing way to think about a 
non~homogeneous Poisson process is to consider the variation 
in the intensity of the event stream over time. The ccncept 
of an intensity function whose integral meets the criteria 
of the definition in paragragh B above is essential tc the 
modeling and simulation of a non-homogeneous Poisson 
process. When one starts with the existence of A(t), the 
A(t) 15 often called the integrated intensity (or rate) 


function. 


The intensity functuion reveals the instantaneous rate 
of arrivals (in the event stream) as a function of time. 
(This function must be a positive function.) For example, 
if telephcne calls arrive at a switchboard at the rate of 5 
per hour at 0900 (9:00 a.m.), increase to a peak rate of 20 
per hour at 1300 (1:00 p.m.) , then decrease to a rate of 5 
Ber hour at 1700 (5:00 p.m.), the intensity function could 
lock something like Figure 2a. Теп, ру рота the 
integral of the intensity function over the interval of 
interest (i.e. from 0900 to 1700) it is obvious that a 
ncnotone-increasing,  right-continuous and bounded function 
is obtained (see Figure 2b). If the assumption is made that 
the arrival strean is a Poisson process, i.e. has 
independent increments, then the number of calls received in 
any chosen interval (e.g. 0900-1000, 1230-1315, 1107-1632) 
1s distributed as a Poisson random variable with parameter 
equal tc the difference bétween the integral evaluated at 
the right end point of the interval and the integral 
evaluated at the left end point of the interval. These 
values may be read directly from Figure 25. Specritrcalle, 


the number of calls received in an eight-hour working day is 


19 


———— ——————— р — — 


a Ecissor random variable with parameter 120. 


Although it 15 unlikely that telephone calls would 
arrive at a switchboard in accordance with the convenient 
parabolic intensity function of Figure 2a, the example 
serves to illustrate two important points. First, it would 
not be realistic to assume that the arrival rate of 
telephone calls at a switchboard would be constant 
throughout a working day. Some function that descrites an 
initially increasing and finally decreasing rate of arrivals 
seems more akin to reality. The importance of being able to 
model a non-homogeneouS Poisson process is therefore 
established. 


Secondly, it is obvious that the definition of a 
ncn-homogeneous Poisson process is not always used as a 
starting point fcr modeling operational processes. Event 
streams are usually thought of in terms of their underlying 
intensity functions. The idea of an intensity function 
applies to any model for an event stream, and is not 
specific to a Pcisson process. The further step of modeling 
the physical process as a non-homogeneous Poisson process by 
assuming that the process has independent increments 15 
taken either on the basis of empirical evidence or physical 
reasoning. Testing for independent increments in a point 
meecess is discussed in COX and LEWIS [ Ref. 3, ch. .6] and 
LEWIS [Ref. 9). The main physical reason for assuming 
independent increments, and therefore a Poisson process, is 
that the cperational process is the superposition of many 
individual event streams. For instance, in a computer 
center equijped with several interactive time-sharing 
terminals, the event stream of users at each specific 
terrinal wight be assumed to be an arbitrary point process 
with a certain intensity function. The event stream seen by 
the central processing unit is then the sum total (or 


Superposition) of the event streams of the individual 
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terminals and, if there are enough terminals, should be 
approximately Poisson. IO opel ty of superposition of 


intensity functions is discussed in the next paragraph. 


A Calls per hour 
20 
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time of day 
Figure 2a - Rate of Arrival of Telephone Calls at a Switchboard 
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Figure 25 - Average Total Calls Received 
Begure 2 - RATS AND INTEGRATED RATS FUNCTIONS FOR TELEPHONE 
CALLS ARRIVING AT A SWITCHBOARD 


D DECOMPOSITION AND SUPERPOSITION OF POISSON PROCESSES. 


To oktain a Poisson process with intensity function 
mete) = a(t) + ào (t); we may superpose two Pcisson 
processes, each of intensity àq (t) and ào (t). 


We may decompose the intensity function of any event 
stream intc twc or more component event streams. However if 
we then superpose the component streams, we will recover the 
original tyre process only if we started with a Pcisson 
process. For example begin with a renewal process with 


intensity y(t) = and let v A If two renewal 


2 * 


processes with intensities v, and v are  superposed, the 


i 2 
resulting prccess is not a renewal process. 


This urique property may be exploited when simulating 
Poisson processes. It permits the partitioning cf the 
intensity function to take advantage of any special 
properties of its component parts. This is the basis for 
the methcd used in the Poisson-Decomposition and Gap 


Statistic Algorithm discussed in latər sections. 


me 2WO EASIC METHODS СЕ GENZRATING A NON-HOMOGENEOUS 
merssON EBOCESS 


ER ЕШ ШЕР ee ee E 2 ee cee ee ee O O NÓ oe ee сш A |. сир Нр > 


Ccnsider the non-homogeneous Poisson process with 
intensity function (Т) > О, 0<t<tor on the interval 


(0,5,1. The integral of the intensity function is then 
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A(t) and the number of events in the interval is a Pcisson 


randon variable with parameter A (t9) - A(0) = m (or 
= 1 = 0 = 1 ж ж 
A (tg) ug Since A (0) fo nit) dit 0). Now if Тї, сез; De 


are events in a unit homogeneous Poisson process (i.e. a 
Poisson process with constant intensity function 
КОМ) = A" = 1) then A (тж), ..- АТТ (тө) are events in the 


non~homogenecus Poisson process. 


This result gives a procedure for simulating a 
ncn-homogeneous Poisson process, starting with a homogeneous 
Poisson process, which is analogous to the probability 
integral transform method of producing random samples from a 
continuous distribution with distribution function F(x) when 
starting with uniformly distributed random variables. The 
latter method is essentially that if the inverse of F(x) can 
be found, then by generating uniform (0,1) variates Use 
E. u the values 


x, = Fu), к. en E (a ) 


ccmprise a sample of variates from the desired distribution. 


The right-continuous monotone increasing function 
A (t) describing the ncn-homogeneous Poisson process on the 
interval (0,tg] can be thought of as a distribution function 
ИАС has been "scaled" by the тас согы (oe nice 
aus.) = їй then ПСЕ - 1, thus AMAL, 15 а valid 
Me tribution function on Ve 


Tc implement the time-scale transformation procedure 
one can use the following basic result for homogeneous 
Pcisson processes: given that n events have occurred in a 
homogeneous Poisson process over a fixed interval  (C,tg], 
those events are uniformly distributed on the interval. A 
Pecot Of this property is given in PARZEN [Ref. 14, p. 190). 


Therefore, if events are generated as a unit Poisson process 
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on an interval of length Up by [iret obtaining a 
Poisson (Hg) random variate n, and then letting the order 
statistics frcm a random sample of uniform (0,10) variates 
be the times to events in the unit Poisson process, and the 
times to these events are then transformed by the inverse of 
the integrated intensity function, i.e. А (о), the effect 
is the same as that obtained from the probability integral 


transforn method. Figure 3 illustrates this method 
graphically and also provides a flow chart. Note the 
difference however between this procedure and the 


probability integral transform procedure. In generating a 
unit Poisscn process by this method we need a sample whose 
size is randcm (and could be zero) i.e. Poisson (Ho) - The 
prcbability integral transform method simply transforms a 
fixed number cf uniform (0,1) variates into variates from 


some other distribution. 


In Appendix A it is proved that scaling both the 
EStripution function F(x) and the interval (0,1) ty the 
Same factcr dces not affect the validity of the probability 


integral transform method. 


The unit homogeneous Poisson process may also be 
generated Ey adding a sequence of unit expenential variates 


(variates frcm an exponentially distributed random variable 


with mean = 1) until the sequence of partial sums of the 
random variatles first exceeds ug (see Appendix р). This 
acccmplishes two things. First it provides an ordered 


sequence cf events from a homogeneous Poisson process on the 
interval (2,101. Second, it determines a realization of the 
Poisson randcm variable Ne, with parameter ÁA(tg) -Ug. Note 
that given n, the times to events are uniform order 


Statistics. 
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The difference between these methods for generating 
a homogeneous Poisson process should also be noted. The 
first method requires a Poisson variate and ordering of that 
humber of uniform random variables. The second requires 
generaticn of independent exponential variates. The second 
method is probably most basic in that it requires only 
exponentially distributed random variables, and these can be 
obtained from uniform variates by the inverse probability 
integral transform, which is just a logarithm. The method 


is, however, not always the most efficient. 





Du Senna ana Order Statistics Algorcitán for a 
Moisson Process 
TAIS method requires the result of a theoren 


sketched by LEWIS and COX (Ref. 3, ch. 2] and restated here 
for convenience and continuity; it is an extension of the 
result on conditioning in a homogeneous Poisson process 
Echo results in a conditional uniform distribution of the 


mmes to events. 


Theorem 1: Assume that a non-homogeneous Poisson process is 


observed for a fixed time (0,tg9], so that the number of 
events in  (0,tg], Neo has a Poisson distribution with 
parameter A (tg) = MD) = Hg 
times-to-events for the process in (0,541, and if Че. = П, 
conditional on having observed n(>0) events in (0,541. the 


SA DO, denote 
Је n 


ENDE are dsstributed as order statistics from a sample with 


As tribution function 


EE E) 
же O) 
nis theorem reduces the problem of simulating a 
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non-homogeneous Poisson process to that cf generating а 
Pcisson number of order statistics from a fixed distribution 
function. That is, given an intensity function over an 
interval  (a,b], whose definite integral A(t) = ГЕ) (в) 46 16 
bounded and right-continuous on the interval, an ordered 


sample, frcm a population with distribution function 


DCUM Ua) 


BEDS). AUS u 


will yield the desired non-homogeneous Poisson process 
defined су the intensity function NE) Quita DOE 
Simplicity the interval will hereafter be assumed to have 
its left end point at zero (a = 0) and its right end point 
at some arbitrary, but fixed point tye (b = 54). Using this 
(9,5, Ј interval results in no loss of generality, and (1) 


becomes identical with the expression in the theoren. 


Many methods exist for obtaining the necessary order 
statistics. Тһе inverse integral transform explained above, 
Mecemposition of the density function (see LEWIS Ref. 8), 
or the rejection-acceptance method discussed later In 


Section IV-D, are all possibilities. 


For the family of intensity functions addressed in 
this thesis the non-hcmogeneous Poisson process is obtained 
by a combination of Poisson decomposition, an algorithm of 
LEWIS and SEEDLER (Ref. 13] for obtaining a non-homogeneous 
Pcisson process with log-linear intensity function, and the 
conditioning and order statistics theorem given above. The 


procedure involves four basic steps. 


1. The intensity function is decomposed into two 


ccmpenents; i.e. Alt) = A(t) + A*(t). 
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2. A gap statistics algorithm is used to generate 
a ncn-homogeneous Poisson process from one of the 
components, A(t), which is chosen such that it has 


a special structure (log-linear). 


See Gegecemon—-acceptance routine 15 used to 
produce a sanple from the remaining component, 
1*(t), i.e. a Poisson number of ordered variates 
are generated. This algorithm is described in 
detail cin SeenON 11 =. This sample, when 


ordered, becomes a Poisson process also. 


4. The Poisson processes of the compcnent 
intensity functions are then superposed to produce 


the desired non-homogeneous Poisson process. 


Figures 4, 5 and 6 (Section III) illustrate the four general 


steps of this procedure. 


Ncte: Theorem 1 provides a second way to generate 
the unit homogeneous Poisson process required by the 
time-scale transformation algorithm previously described. 


First we may generate a variate from the Poisscn random 


variable Meg with parameter но, say Mea Ds Then 
conditional upon М. = п, п uniform order Statistics can be 
generated cn the interval. For reasons explained in 


Appendix D, this second method was used in the computer 
program implementation of the time-scale transformation 
method. 
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ПИПИНА ОНЧНАТЕЕОН So LECTION OF DEGREE-TWO EXPONENTIAL 
POLYNOMIAL INTENSITY FUNCTION 


The intensity function identifying the non-homogeneous 
Poisson process investigated in this paper is of the form, 
62) 


А СЕ = exp (ag det P 


l 2 


where а a, and a. are real constants. 
0 1 2 

This intensity function was selected for three reasons. 

First, an intensity function must always be positive (or 

zero) if it 1s to be meaningful. The above function is 


0 Ш 
nas intensity function, by proper choice of constants, can 


non-negative for all values of а, A and y Secondly, 
be used tc represent the four different types of event 
streams mentioned previously in Section I. And finally, the 
selection or this intensity function leads to Simple 
statistical procedures; Monde cas, “сес LEWIS Бег, 9, 


2-30—34). 
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A. GENERAL 


Assuming an intensity function of the form 


A(t) = exp(a, + a, + at?) 
0 1 
three algorithms ТОП generating the corresponding 
ncn-homogeneous Poisson process are discussed. These 


algorithms are based on the two general methods presented in 
Section II-E, and the decomposition and superposition 


property of Poisson processes. 


Eee CEME=-SCALE TRANSFORMATION ALGORITHM (ALGORITHM А) 


1. Step One 


By definition of a non-homogeneous Poisson process, 
the total number of events observed over a fixed interval 
(0,50) is itself a cc distributed random variable, Nes 
with parameter Mo = 497 A(t) at. The first step cf ave 
algorithm is to determine the value of the parameter Ug · 
Although an explicit, closed-form expression for the above 
integral cannot. be found, a series representation does 


exist. Except POL a constant tacto. this series 
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representaticn assumes the form cf the error functicn cr of 
Dawson's integral. A negative value for the coefficient of 


the second degree term, a, yields the error function forn: 


a^ ty 2 
-u 
ug 7 x ( 27 2 u: a du) 
Ут 0 Ут 0 
whereas a positive value for results in the Dawson's 


integral fors: 


T 
O 
|| 
N 
ct 
N 
КОЛ ГУ 
ct 
N | 
N 
0 
N 
(D 
с 
~ 
as 
С 
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ct 
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ct 
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In the above expressions K, t and € are uniquely 


Es 2 


determined by the coefficients a and the end points 


(ua 
of the interval over which the intensity function applies. 
A detailed derivation cf above relationships is given in 


Appendix C. 


Evaluation of the error function 25 а FORTRAN 
supplied procedure and requires only that the proper 
arguments be calculated and provided to the FORTRAN FUNCTION 
ЕНЕ or CERF [Ref. 15]. Evaluation of Dawscn's integral is 
best accomplished through use of the IMSL (International 
Mathematical and Statistical Libraries, Inc.) FUNCTION MMDAW 
(Ref. 6]. The accuracy of the function values calculated by 
these Feutines is linited only by the precision 


emaracteristics of the computer. 


2. Step Two 


om ou рар 


Cnce the parameter of the Poisson random variable 


N is determined, a realization on that ranđom variable is 
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required. (The approach is somewhat backwards since it 
first determines hcw many events occurred over the interval. 
It then distributes that fixed number of events over the 
interval in accordance with the non-homogeneous Pcisson 
process described by the intensity function. The importance 
of Theorem 1 now becomes evident since it assures the 
validity of such a procedure.) Generation of Poisson 
variates, especially those with large parameter values, is a 
Бешріех procedure in itself if efficiency in terms of 


computer time and memory requirements is desired. This 


prcblem is discussed later in Section IV-B. For the 
present, assume that the requisite variate has been 
prcduced, i.e. N = п. 

EQ 


3. Step Th 


Given that n events have occurred over the interval 
(0,6 ]) we then distribute п events along an interval of 
length Ho in accordance with a homogeneous Poisson Prccess. 
Since events іп а homogeneous Pcisson process are uniformly 
distributed over an interval (given that n events have 
occurred), this step merely requires that n uniform (0,1) 
variates be generated, ordered.from lowest to highest, and 


then each multiplied by the factor u The values in this 


Qs 
n-element vector, (ul, чл, Lu ut): correspond to tie 


Euues plotted on the vertical axis in Figure 3. 


Each event in the homogeneous Poisson process must 
be transformed by the inverse of the integral of the 
intensity function. Letting gt A (s)ds = A(t), the inverse 
À 7i (e) applied to each event E the homogeneous Poisson 


process, will produce a corresponding event in the 
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non-homogeneous Poisscn process, Jes АЁ (а 
-1 
A (о.) 


integral of this specific intensity function cannot usually 


ee 
ZU ес. The difficulty is that since the 


be explicitly expressed, the form of its inverse usually 
eludes any convenient computational formula expression. The 
unigue pesiticn on (29,521 the inverse determines for each 
input value can be found to any degree of accuracy desired, 
by iterative, numerical methods. The Newton-Raphson method 
is easily employed and very efficient in the present 
scenario. Its implementation is explained in Section IV-C. 
Since the function A(t) is strictly monotone increasing, 
the inverse function A71 (u) applied to an ordered sequence 
of input values results in an ordered sequence of output 
values. Therefore, the toe wE у = are the times cf events 
in the non-hcmogeneous Poisson process and the algorithm is 


complete. 


C.  TIME-SCALE TRANSFORMATION ALGCRITHM, ALTERNATE 
(ALGORITHM A!) 


An alternative approach to the time-scale transformation 
method described akove 1s to generate the required 
hcmcgeneous Poisson process by using the fact that in this 
process the random times between events are independently 
exponentially distributed. Thus one generates unit 
exponential variates until their sum exceeds М0” Тһе 
partial sums give the times to events and the number of 
partial sums less than or equal to uo is a Poisson (ug) 
variate. Note that the Poisson variate comes out as a 
by-product in this procedure rather than as a pre-product as 


in Step Two of Algorithm A above. 


Althcugh this methcd combines Step Two and Step Three of 


Algorithm A into a single procedure, it is not necessarily 
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the best method. Because it requires the use of the 
additive method of Poisson variate generaticn, it becomes 
inefficient for Poisson processes with many events (see 
Appendix D). For this reason Algorithm A instead of 
Algorithm A* was used when implementing the time-scale 


transformaticn method into a computer program. 


р.  POISSCN-DECOMPOSITION AND GAP STATISTIC ALGORITHM 
(ALGORITHM B) 


It is recommended that the reader refer to Figures 4, 5 


and 6 for a better understanding of the following steps.) 


The Poisson-decomposition and gap statistic 
algorithm begins with an examination of the coefficients of 
the intensity function. By doing so, the intensity function 
meme categorized into one of six possible configurations. 
These six cases are discussed in LEWIS and SREDLER 
ER sr. 11]. Examples of each case are illustrated in Figures 


 спголаһ 12 located at the end of this section. 


2. Step Iwo 


а ОТ cher intensity function A(t) is monotone 
increasing or monotone descreasing over the interval (Cases 
MEL, IV and V; see Figures 7, 8, 10 and 11) the intensity 
function is decomposed into two separate intensity functions 
Over the same interval. The resulting intensity functions 


РЕС of the £cram; 


| >= 
сї 
w 


= exp (Y y ш Ye 
and 
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t +a t^) - exp lyy + Y,t) 


AF(t) = A(t) = Alt) = expla, + а, 2 


iis Gleammethat A(t) + A*(t) = A(t). 


De If either of the two cases not covered by 2a 
abcve occurs, à(t) will be monotone increasing (decreasing) 
on the subinterval (0,b] and monotone decreasing 
(increasing) on the complementary  subinterval  (b,tg] (see 
Figures 3 and 6), where b 15 а unique point within the 
interval at which A(t) has a maximum (minimum) value. By 
dividing the interval properly into two disjoint, contiguous 
subintervals, each subinterval may be treated as explained 
mn 2a. Subsequent steps are applied to each of the two 


intervals separately and the results combined. 


An Shere ient algorithn for generating a 
ncn-homogeneous Poisscn process with a log-linear intensity 
Mugection (1.6. A(t) = exp (8, + 3,4)) is presented by LEWIS 
and SIEDLER fRef, 13). This algorithm generates the 


non-homogenecus Poisson process through the use of gap 


Bea tistics. (A comparison of the gap statistic technigue 
with the conventional integral transform technique is 
discussed in Appendix C.) By judicious selection cf the 


coefficients cf the log-linear intensity function, most of 
the total area under the original intensity function 4) (t) 
will be contained under the function A(t). Therefore most 
of the events in the non-homogeneous Poisson process with 
intensity function A(t) can be accounted for by employing 
the gap statistics algorithm on the intensity function 
pt). 
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Е I NN Ea 


A us 


0 1 2 


Step 3 - Generate events T; from log-linear intensity 
function, A(t), using gap statistic algorithm. 
Events will be ordered. 
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A 


0 1 2 

Step 4 - Generate events 5. from intensity function 
A*(t) = A(t) - A(t) using rejection algorithm. 
Events will not be ordered. 


0 1 2 
Step 5a - Order events from Step 4, 


Шы 7 lg "s Но 
0 | 2 


Step 5b - Merge (superpose) events from Steps 3 and 5a, 
Result is event stream Ty; 1, МУ РОТ 
ЕЕ) (BE 


ПЕ от ОТАСДАЙ (CONTINUED) 


з 
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All that remains to be done is to generate an 
ordered sample of some given size, on the interval (0,5,1, 
frcm the remaining component of the original intensity 


Bunction, i.e.  x*(t). 


The number of events in the interval is a Poisson 
random variable КЕ 7 with parameter ut = s 001% (t) dt. 
Once a realization on this random variable is obtained, i.e. 
Ni =n', Theorem 1 may be invoked. Since A* (t)/ u9 is 
more easily evaluated than its indefinite integral, the 
rejection technique (explained in Section IV-D) is used to 
generate the n* required variates. The rejection technique 
is not, in general, always an efficient method for variate 
generation unless great care is taken. Yet in this 
deccmpositicn scenario, it will be used to generate cnly a 
small percent of the total events required by the original 
intensity function AENA The Majority of the events will 
be generated by the efficient gap statistics algorithm. The 
efficiency gains should more than compensate for any 


efficiency lcsses due to the use of the rejection technique. 
5. step Five 


The events produced by Step 3 will be in order on 
me псехтај (9,60). The events produced in Step 4 will not 
be in order cn (0,tg]. By ordering the events from Step 4, 
(which are few in number compared to the total number of 
events in the non-homogeneous Poisson process) it is 
possible tc superpose the two ordered event streams. The 
merced event streams produce a new event stream from the 


original intensity function Roce es 
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TY.: ALGORITHM IMPLEMENTATION 


A A A O MA E O O о ee 


A. GENERAL 


This section explains the several specific techniques 
that were applied to implement the two competing algorithms 
(MWigoritha. A and Algorithm В) into PORTRAN computer 
pregrams. Detailed discussion of the various subprograms is 
avoided since References 11 and 13 and attached program 
listings provide such information. The Algorithms A and B 
have some subprogram requirements in common while other 
subprograms are unique to one algorithm or the other. This 
section will discuss first those requirements common to both 
algoritho<, then those needed only by the time-scale 
transformation algorithm (Algorithm A) and finally those 
unique to the Poisson-decomposition and gap statistic 
algorithm (Algorithm 8). Hereafter differentiation between 
the algorithms and the FORTRAN computer progran 
lmplementaticn of the algorithms will not always be made. 


The meaning will be clear írom the context. 


See COUMCN REQUIREMENTS 


1. Integration of a Deqree-Two Exponential Pclyromial 


Í -a e ес — ap SSS эшш — ee сс с. 


Pumction cver a Fixed Interval 


о => ОНЕ se SSE сар A шшш > zum eb ri mpm O AO E шшш O. ee ee b 


Algorithm A requires that the intensity function 
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А (+) be integrated cver a fixed interval in order to 
determine the value of the parameter of the Poisson random 
variable that governs the number of events that are observed 


to occur in the non-homogeneous Poisson process. Algorithm B 


MESES tide the intensity function aA*(t) = lt) = Alt), 
be integrated over a fixed interval for the same reason. 
Since 

b b b 

ME ОЕ = еа у ағ 

а а а 


and A(t) has an explicit expression for its indefinite 
Emtegral, i.e. 


b exp (Y4t) 
ME (ја = = exp (Y g) — 
a а 


the problem fcr both algorithms is reduced to computing the 
value for 
le b 2 


ete cee | exp(a, + a,t + a,t”)dt 
a a 


EHEROUTINE HELP employs IMSL FUNCTION MMDAW or the FCRTRAN 
supplied procedure DERF, as appropriate, to perform this 
Calculation. Section III-B and Appendix B discuss how this 
computation could also be made using a ccnvergent series 


representaticn. 


2. Generation of a oisson Variate with a Given 


Parameter 


Both candidate algorithms require at least one 
meer zaticn on a Poisson distributed random variable with a 


given parameter. The value assumed by the random variable 
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is the number of events observed in a specific interval of 
time over which the occurrence of events is governed by a 
given intensity function (і.е. A(t) for Algorithm А and 
A *(t) fcr Algorithm B). The nature of most real event 
streams which lend themselves to analysis using a 
non-homogeneous Poisson process model is that they consist 
of a large tctal number of events. This will be the case 
either if a dense process is observed over an interval of 
shcrt or mcderate length or if a "sparse" process is 
observed over an extended interval (see the arrivals at an 
intensive care unit given in LEWIS (Ref. 9]). The point is 
that both algorithms must be flexible enough to generate 
ncn-homogeneous Poisscn processes that result in high 
numbers cf events occurring. Since the number of events 
occurring over any interval in such a process is a Poisson 


random variable, it becomes necessary to be able to generate 


a variate from a Poisson distribution with a large 
parameter, i.e. large mean. The two most theoretically 
straightfcrward algorithms for generating Pcisson 
distributed variates (the additive and multiplicative 


methods) Eeccme computationally cumbersome as the parameter 
of the random variable increases. Both the additive and 
multiplicative algorithms and the practical deficiencies of 


each are discussed in Appendix D. 


The technique selected to deliver a Poisson variate 
On demand tc both Algorithm A and Algorithem B is the Gamma 
method. It is developed and explained in an unpublished 
book by AHRENS and DIETER (Ref. 1]. A paraphrased account 
of the development of this algorithm is included in 
Appendix Г. The main advantage of AHRENS and DIETER's Gamma 
method is that whereas the computer time required for the 
additive and multiplicative methods is proportional to the 
parameter of the Poisson distribution being sampled, the 
computer time required by the Gamma method is propcrtional 
to the lcgarithn of the parameter. The SUBROUTINE PVAR 
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employs the Gamma algorithm to return a Poisson variate when 


given a parameter as an input. 


3. Event Storage 


Any algorithm which simulates а non-homogeneous 
Poisson process must have some mechanism to provide the user 
with information that completely describes the 
realization(s) of the process simulated. Since the specific 
arrangement of the stream of events over the interval is the 
informaticn cf interest to the analyst, the locaticn, or 
time of occurrence, of each single event on the interval 
must be stored. Equivalently, the spacings between events 
would ccmpletely describe the realizations on a 
non-homogeneous Poisson process. Thus event spacing 
informaticn rather than event location information could be 
Stered. (Programs written for candidate algorithms A and B 
both provide the cption for the user tc demand either event 
location or event spacing information.) When using either 
algorithm, an array large enough to hold location or spacing 
data for each event generated by the algorithm must be 
created. Since the number of events observed in any 
realization of a non-homogeneous Poisson process is itself a 
random variable, the precise size of the array cannot be 
determined a priori. If the programs are to have any value 
for general application, they must be able to accept 
intensity function parameters which will demand large 
numbers of events when simulated. Using the somewhat 
arbitrary assumption that an event stream with an average of 
4500 events is sufficient for most Simulation scenarics, a 
fixed array with a capacity for 5000 events is used in the 
pregrams implementing both of the algorithms. If the value 
of the intensity function integrated over the interval is as 
high as the maximum limit of 4500 (i.e. the number of events 


in the interval is Poisson distributed with msan = 4500) 
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then the array of length 5000 allows the Poisson random 
variable to exceed its mean by 7.45 standard deviations 
before an array overflow is encountered. This ау 
unlikely event is of such rare occurrence (less than one 
chance in a billion) that it may be disregarded. However, 
should it cccur, the programs will reinitialize themselves 
and generate a new Poisson process. Also an error indicator 
is incremented. This error indicator (IER) may be written 
on demand. Its value is the number of times the program was 
forced tc akort generating a Poisson process, reinitialize 


itself and start again. 


A 5000 element array is small enough to avoid its 
being an undue memory requirement burden on most operating 
systems, yet large enough Lo accommodate nost 
non-homogenecus Poisson processes of interest. Its choice, 
though arbitrary, was based upon the above two 


considerations. 


SPECIAL REQUIREMENTS OF THE TINE-SCALE TRANSFOREATICN 
ALGCRITHMÁ (ALGORITHM A) 


As explained in Section III-B, once the number of 
events ckserved in the non-homogeneous Poisson process is 
established (i.e. “<n ПЕ ts necessary to construct" a 
homogeneous Poisson process consisting of n events over an 
interval cf length ug units. This is easily done by 
generating n uniformly distributed variates on the interval 
(071) and then scaling each by the factor Hg * The LLRANDOM 
computer likrary package developed by LEWIS and LEARMONTH 


[Ref. 12] is a very efficient source of such variates. Once 
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an array of n uniform (0,1) variates is obtained from 
LLRANDOM, each element of the array must be multiplied by 
the apprepriate scaling factor. It is these resulting 
numbers which must be acted upon by the inverse of the 
шеста са “intensity function to yield the location of 
events in the non-homogeneous Poisson process on the 


original interval (0,t,]- 


2. Sorting of Events 


Tc te of much practical value a Simulation routine 
for a non-hcmogeneous Poisson process should order the 
events in the interval from first to last. In Algorithm A, 
this ordering may be done before applying the inverse 
integrated intensity function to the elements in the 
homogeneous Poisson process. The monotonic nature of the 
integrated intensity function maintains the relative order 
of all elements after they have been transformed by the 
inverse function (see Figure 3). Of course the ordering 
could be dcne after the transformations also. Hem 
implementing Algorithm A, the uniform variates were scaled 
by the factor ug 


before being transformed by the inverse of the integrated 


then ordered, from lowest to highest, 


Entensity function. 


Ordering of large arrays of numbers 15 а time 
consuming operation on the computer. There are many 
ordering algorithms of varying degrees of sophistication and 
efficiency. Because ordering is unavoidable when using 
Algorithm A, selection of an efficient ordering routine is 
ENortant. The ordering algorithm used in this 
izplementation was the W. R. CHURCH computer center library 
routine FXSORT which employs Singleton's version of the 
partition exchange sort. (program listing of PXSERT Is 


provided in the computer center's catalogue of library 
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routines.) The PXSORT routine appears to be the most 
efficient crdering routine readily available Tor the 
purpose. (It is acknowledged that more efficient routines 
specifically tailored to the problem of ordering uniform 
random variates may possibly have been developed that would 
improve the overall efficiency of the program implementing 
Algorithm A.) 


CERO o CER. AA | Í T с җы 





The essence of Algorithm A is the application of the 
inverse cf the integrated intensity function to each event 
in the homogeneous Poisson process on the interval (C,ug]- 
AS menticned before, the intensity function, A(t), under 
consideration does not yield an explicit expression fer the 
integrated intensity function, A(t) = F(t), that must be 
inverted. (Note: F(t) is a "scaled" distribution function). 
Hence neither does a ccmputationally convenient expression 
for this inverse function exist. Numerical methods must be 
employed to transform the position of each event on the 
interval (0,1, ] in the homogeneous process tc its 
corresponding position on the interval (0,t,] over which the 


Simulated non-homogeneous Poisson process is to be produced. 


The Newton-Raphson technique was used to accomplish 
the transformation. This technique allows for iterative 
approximations which converge to the true transformed 
a where En = етеу The 
iterations continue for each value u, until an estimate t' 
Ber its corresponding EE is found that satisfies 
ФЕ”) – о. | < = , where © is a predetermined tclerance 


тре == 
level. (The specific value for © used was Uy x 10 .) 


values. (i.e. tie tor e e.p Е 


In the Newton-Raphson technique, given a function 


ES tne objective is to find the solution, x*, satisfying 





the equation h(x*) = 0. Letting xy be the initial estimate 
for x* (based upon some reasonable criteria) the next 
estimate for x*, that is, х, can be calculated from the 


fundamental iteration relationship, 


h(x,) 


"EL ATK i 


Kell 


This assumes that the function 1(•) and its derivative h'(e) 
can be evaluated at xy. The process is repeated k times 
Peter |5 (х1) 1 <E ; or, until ІХхұ - хұн! Ge 2 This 
procedure is nothing more than using the first two terms of 
the Taylcr series expansion of the function h(e) to locate 
the x-intercept of the line tangent to h(x) at the point 
h(x,). That intercept value becomes the next approximation 
for the roct of h(x) and the procedure is repeated. A 
graphical explanaticn of the Newton-Raphson method is 


presented in Appendix E. 


Applying the Newton-Raphson method to the present 
problem requires some special modifications. The problem is 
to find a value ti such that P (ti) =u, where u, is known, 
е. to Im DLE (u,). Direct application of the 
inverse F (») is impossible since its functional forn 


defies all but the most abstract mathematical expression. 


However if a new function G; (t) = F(t) - u, is defined and 
if its root t*, determined (i.e. a t* such that Gi (t*) = 9) 
then F(t*) - ab = 0, and F(t*) = а... TRUS ue -1 is the 


value desired. 


In order to apply the Newton-Raphson method to the 
Bunction G; (t) it must be possible to calculate both G; (t) 
and its derivative. Since G. (t) di ferc rron F(t) only by a 
constant, the problem reduces to that of evaluating F(t). 
This has already been done, as previously explained in 


Section III-B, and merely requires the assistance of 
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SUEROUTINE HELP. Now G * (t) = Р*(+), so taking the 
derivative cf G, (t) returns the intensity function  A(t) 
which is easily evaluated for any t. Clearly, the 


Newton-Raphson method may be used. 


It shculd be noted here that for each u, there is a 
correspondinc function G , (t) cn which the Newton-Raphson 
technique must be used. Since each use of the 
Newton-Raphson method may require several iteraticns to 
arrive at a value t* which satisfies the tolerance 
criterion, Pt is readily evident that for large a, 
considerable computational effort is reguired to obtain all 
the values ti MM 1^, M Thesmumber of iterations 
required to arrive at a suitable approximaticn for each ti 
15 hignly sensitive to the accuracy of the initial 
approximation (designated E If the initial approximation 
is close to the actual tir the Newton-Raphson method will 
converge very quickly. тї де ea rra pproximati1cn".1s 
pocr, ccnvergence could be much slower. The procedure used 
to select initial approximations for each ti will therefore 
have a profound effect upon the overall efficiency of 
Algorithm A. 


Selection cf these initial approximations is done by 
partitioning the interval (9,%,1 into equal length segments. 
The number of segments is equal to min(10, n/4]. The 
function F(e) is evaluated at the end points of each segment 
and these end points, with their corresponding function 
values are stored in an array. This procedure is performed 
O SUBROUTINE  BNCHMK. See Appendix E for a graphical 


representaticn. 


Ша ап ап села! spproxismatron for the t; which 
corresponds to any given Use the array of function values is 
searched until two adjacent function values which bracket u, 


i 
are located. These adjacent function values uniquely 
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identify the segment of the interval 1 in which the 
elusi ve ti is located. Either end point cf the segment 
would serve as a good initial approximation t; „for the 
Newton-Raphson method. However, for the purpose of this 
implementaticn, the end point which yielded the function 


value clcser to u; is used for the initial approximaticn. 


The decision to divide the interval (0,tg] into 
шіп (10, n/8] segments was based upon empirical results. Of 
Several proposed segmenting schemes that one which resulted 
in the fewest total number of distribution zunetıon 
evaluations over several intensity function/interval 
Scenarios was chosen. Higher values than n/4 for sparse 
event streams may produce better results. But n/4 gave good 
results for dense streams and adequate results for sparse 
streams. It appeared to be the best ccmpromise asa 
candidate fcr general usage. (The number of function 
evaluations of Gi(t) for each t; averaged between 2.2 and 


У тог this partitioning scheme.) 


Alternative methods would be to use for the initial 

approximation for t one of the following values: 

- 4-1 (<t i) 

" pet on 
Neither cf these other methods were attempted so it is 
uncertain whether they would be more or less efficient than 
the method which was used. Efriciency differences among the 
three methods are probably not substantial since each will 


usually yield a good first approximation. 
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ЕИО СЕЕ REQUIREMENTS OP THE DECOMPOSITION ALGORITHM 
(ALGORITHM B) 


I Intensity Function Categorization 


=F ee coe ш жне жы шт A сым 


The decomposition routine selected depends on which 
of six possible shape categories the intensity function 
À (t) falls into (cor. Figures 7 ~- tirough: 1276 An 
examination cf the constants 20, Ay and а. in the intensity 
Function and the interval (0,tg] over which the 
non-homogeneous Poisson process is to occur will uniquely 
identify the category of the intensity function. A thorough 
discussicn of this procedure is presented in Ref. 11, and is 
nct rerroduced here. Implementation of the category 
identificaticn procedure requires a lengthy sequence of 


decision statements within the computer progran. 


; Selection of the Imbedded Log-Linear Intensity 


eS eee хо 





Once the intensity function has been categorized at 
must be decomposed in accordance with a complex scheme 
described in Ref. 11. The objective of the decomposition 
scheme is te fit a log-linear curve (or curves) completely 
underneath A(t) on the interval (i.e. E SO 0 <tSto) 
in such a way as tO maximize the area under the log-linear 
curve(s). This 1s done by partitioning of the interval 
(0,t o] Moto sub’intervals (0,5) and (b,tg] if necessary, and 
then by ctroper selection of the coefficients Yo and ҮТ TOT 
the log-linear ето (5) These coefficients are 


functions cf the coefficients Uy” а. and e ia the original 
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intensity function A(t) and of the interval (0, tg ]. As the 
program advances through the categorizing decision 
statements, the proper coefficients are computed fcr the 


appropriate Iicg-linear function(s). 





The precursor to LEWIS and SHEDLER [Ref. 11] was a 
paper by the same authors (Ref. 13] proposing a gap 
statistics algorithm for simulating a non-homogeneous 
Poisson process Wierda og-isnear intensity function 
\ (5) = exp(8g % 815). (Reference 11 reviews this algorithn 
in detail.) As previously noted, Algorithm B divides the 
intensity function AA COSA sum Of two intensity 
munctjons, À(t) and СЕЈ The new intensity function 
Is Chosen to take on the form of a log-linear function 
so that the cap statistics algorithm may be used to generate 
a stream of events from this portion of the original 
intensity function A(t). The SUBROUTINE NHPP2 implements 
E LEWIS and SHEDLER gap statistics algorithm for the input 
Mi tensity function A(t) and returns an appropriate event 
stream to the calling program. Since the use of the gap 
statistics algorithm within Algorithm B is the rationale for 
suspecting it to bea more efficient algorithm than the 
time-scale transformation, it is instructive to examine the 
efficiencies gained in using the gap statistics algcrithn 
Vice a time-scale transformation algorithm on a log-linear 
intensity function. This question is addressed in Appendix 
ci It was found that use of the gap statistics algorithn 
resulted in a reduction in computer time of approximately 
50% from that required by the time-scale transformation 
method. 
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The дар statistics algorithn has efficiently 
produced an event stream from a non~homogeneous Poisson 
process defined by the log-linear intensity function 
x (t) = EXP {Yọ + 715) · But the process to be simulated has 
intensity function x(t) = exp (ag + aq € + БОЕ): ler is 
therefore necessary to superpose an event stream frcm the 


intensity function 


A*(E) = Alt) - Alt) 


2 
exp(a, + a,t + a) ) - exp (Yg T Ye) 


0 ү 


onto the event stream obtained from the log-linear intensity 


unction. 


Once it is determined how many events are to occur 
over an interval with intensity function AE), Muse: 
М, = пў, it is necessary tc select an ordered sample of n' 


to 
variates from a distribution with density function 


le) 

to 
e 
0 


E 





The value n’ will be small compared to the total number of 
events in the original non-homogeneous Poisson process. 
Therefore the need for realizing efficiency in generating 
these n' crdered variates is not usually as crucial to 
overall algcrithm performance as is the need for efficiency 
in producing the non-homogeneous Poisson precess from the 
MaS пеаг intensity, A(t), in Step 3. However if the 


technigue for obtaining these n' ordered variates is 
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extremely inefficient, much or all of the efficiency gains 
realized frcm Step 3 above could be lost here. (In fact an 
example of such a loss of previously gained efficiency is 


documented in this thesis; see Section VI-B.) 


The rejection methcd is particulary useful for 
generating random variates from populations with continuous 
densities that are bounded and which are concentrated on a 
finite interval; this is the case for f*(t). The rejection 
method and an algorithm for its implementaticn are presented 
in LEWIS and SHEDLER [Ref. 11]. A geometric argument will 
suffice to exhibit the principle involved. Consider the 
density functicn £(x), for а random variable X on the 
interval (a,b) іп Figure 13. The maximum ordinate of this 
density 15 с. If another function g(x) = c* > c is 
constructed then the density function f(x) is enclosed 
Arena the rectangle (a,0), (b,0), (a ,c*), (DE) If a 
pcint within this rectangle is selected at random it will 
fall either under the density function or above it. ШЕ 
is under the density curve the abscissa of that point is 
accepted as a variate from the population. TER Ce ne 
lies above the density curve that point is rejected and 
ancther point within the rectangle is selected at random. 
The procedure continues until п! variates have been 
preduced. Тһе random points within the rectangle are easily 
produced Бу generating two independent uniform (0,1) 
variates and scaling them properly resulting in a point 
KEY), where х = а + u (D - а) and y - 1-с%. Then ir 
f(x) 2 y, x is accepted as a variate, otherwise it is not. 
After n! variates have been obtained, they are ordered to 
yield an event stream for a process with intensity function 
ПЕК). SUERQUTINE REJECT employs this method in the 
puegram for Algorithm В. 
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= * 
B > qc 
Ze in | 
| density f(x) | 
| | | 


^2 = D 


The density f(x) is completely enclosed by the rectangle 


a E AAA a A 


Procedure: 
1. Generate 2 independer“ uniform (0,11 variates а, 


апа 12. 
2. Compute: x = a + Duy Vv a yori 
Bem CM d y). 

4. If the point (x,y) lies under the density f(x) accept x 


as a variate; otherwise go to Step 1. 


Example: In the graphical example above X1 wouid be accepted 


as a variate whereas x, would not be accepted. 


2 
(Note: Ideally, c* =c and a and b correspond 


to the lateral limits of the density “(x). This 
minimizes vejection region.) 


ПИ си NNI OREJECTLION-ACCEPTANCE METHOD ОТ ТАЗТАТЕ 


H 
GENERATION EROM AN ARSITRARY DENSITY 
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ӘНШІСІ Пе probability that a "Point on the abscissa 
will not te accepted as a variate is proportional tc the 
area within the rectangle that is not under the density 
function, it is advantageous to make this area as small as 
possible. Therefore it is best to set c* = c and set the 
values a and b equal to the end points of the interval over 
which the density function occurs. This is desirable but 
not necessary. Figure 13 illustrates the more general case 
where the rectangle is larger than necessary. One may be 
required to select a c* > c if c is difficult to determine. 
Seldcm would a and b not coincide with the lateral limits of 


the density however. 


It is obvious from Figure 13 and the preceding 
discussicn cf the rejection nethod of variate generation 
Це it is the "shape" of the density function which is 
critical to the validity of the method and not the fact that 
it integrates to one. Therefore any scaled density would 
preserve the relative shape of the density function. As 
long as the value c* is at least as great as the maximum 
value of the scaled density, the rejection method will yield 
valid variates. The intensity function A*(t) may be thought 
of as a density scaled by the factor Ho = гж (+)а+. 
Mmegerithm B uses the intensity function A*(t) rather than 
pue density function л (6) Ире Pathe were jection routines, 
This eliminates a division step for every function 


evaluaticn required by the method. 


=e Merging Event Streams 


The event streams from the intensity functions л (+) 
and A*(t) must be merged to produce an event stream from 
y (t) = A(t) + à *(t). In the present case there are two 


event streams (one long and one short) which are already 
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ordered. This structure allows for superposition rather 
than a more general (and more time consuming) ordering 
procedure. To do this merging job SUBROUTINE COLATE is 
called. 


E. SUMMARY 


This secticn has presented a general discussion of the 
more significant and interesting procedures used to 
implement Algorithms A and B efficiently in FORTRAN. The 
Mmeader is referred to LEWIS and SCHEDLER [ Ref. 11] for a 
detailed listing of steps for Algorithm B. Program listings 


may also ke consulted. These may be found after Appendix E. 
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A. GENERAL 


Conclusicns drawn from the results of comparing the 
efficiencies of two computer programs are only as sound as 
the measures of effectiveness upon which they are tased. 
These measures of effectiveness are always somewhat 
arbitrary, so their selection should be based upon 


compelling reasoning. 


Even after reasonable methods of effectiveness have been 
defined, only a finite number of trials (usually a small 
finite number) for a given set of circumstances may bə run. 
Therefore, general statements such as "this algorithn is 
better than that algorithm" can seldom be made with complete 
confidence. Yet if the incomplete information gathered 
reveals certain trends, generalizations hypothesized from 


such trends may warrant a high degree of confidence. 


This section describes how the two algorithms are 


ccapared and the rationale for choosing the method employed. 





MEASURES OF EFFECTIVENESS 


Цио пр лас опа:  =рева 


Perhaps the most popular efficiency measure for 
ccmputer programs is their speed of execution. Speed is a 
natural standard as long as computer time remains a costly, 


оаза гсе Lesource. 


Algorithms A and B were compared in terms cf the 
mean time required to generate event streams from a given 
set of intensity functions. It should be noted that Since 
the number of events in each event stream is а random 
variable, the time reguired to generate one event stream is 


also a random variable. 


Several event streams with different intensity 
functions, each having a different expected number of 
events, are produced. Event streams for each of these 
intensity functions are replicated many times by each 
eegorithm. The total execution time required for all 


replicaticns is the "speed" measure of effectiveness. 


The number of replications varies with each 
intensity function. For example, an event stream with an 
expected number of events less than 200 was replicated 100 
times whereas event streams with over 1000 expected events 
were replicated only 30 times. In both cases the total 
number of events produced was of the same order OE 
Magnitude. A priori it seemed reasonable to assume that 
computer time required to execute each algorithm's  progran 


would be roughly proportional to the total number of events 
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generated. Therefore in designing the experiment the 
product cf the expected number of events times the number of 
rerlicaticns (EE N eo ] x * reps) was kept roughly constant in 
order to keep demands on computer time reasonably under 
control. (Note: Program execution times were isolated from 
compilaticn and linkage times for the purpose of the 


computational speed comparison.) 


2. Computer Memory Requirements " 





A seccnd natural efficiency measure for computer 
ргсагапв 15 their core memory requirement. Both programs 
were written with the goal in mind of conserving as much 
memory as possible without unduly increasing program 
complexity. Core requirements reductions could most likely 
be accomriished in both programs through more sophisticated 
programming methods. Therefore this comparison has a major 
weakness as a measure of effectiveness. Recognizing this 


caveat, memory requirements for each program were compared. 


Be Fidelity 


The question of paramount importance 15: How well 
does the program simulate the true process? This is a 
difficult question to answer when dealing with a finite 
number cf realizations on a process whose theoretical 
properties can be validated only after an infinite number of 


realizaticns. 


Several methods of comparison are possible. The 
Simplest is that cf comparing the sample means and sample 
variances of the random variable whose realizations are the 
number of events generated on an interval by а given 


intensity function. This random variable should be  Pcisson 
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distributed with mean equal to the integral of the intensity 
function over the interval. For a Poisson random variable, 
the variance equals the mean. This test, although useful 
for determining if either algorithm varies grossly from 
expected performance, yields no information about hcw the 
algorithms are distributing their simulated events over the 
interval. It is necessary to look at discrete segments of 
the interval and examine the distribution of the number of 
events produced in each segment by our simulated 


non-homocenecus Poisson process. 


From the definition of a non-homogeneous Poisson 
process (see Section II-B) we know that the number of events 
observed over any segment of the interval must be a Poisson 
random variable with parameter equal to the integral of the 
intensity function evaluated over the interval. By dividing 
into several segments the interval over which the Pcisson 
process is being simulated, the number of events in each 
segment can ke observed. Two hypothesis tests may then be 


made for each segment. 


Ше каа = test Ше the dispersion test [Ref. 3, 
oO. 1581. This test seeks to ascertain if the number of 
events observed over each segment is distributed as a 
1' 2! EE n, be k 
observations on the number N of events occurring in a given 


Poisson random variable. Let n n 
| ре 
fixed segment р, and let RE: (n, * +. .+n,)/K. Ir UT ls a 
Poisson distributed random variable, then the statistic 
k - 2 
а 607590) 
qu 7 р, р 
n 
p 
has a distribution which is approximated by a chi-squared 
distribution with k - 1 degrees of freedom. If the interval 
has been partitioned into m segments, then by simulating the 


ncn-homogeneous Poisson process k times, we can perform m 





hypothesis tests at a selected significance levela. The 
test statistic 35 would then be compared to the (1-а) 
percentile of the chi-squared distribution with k- 1 
degrees cf freedom. For a large number of hypothesis tests, 
say 1000, we would expect approximately 1000q rejections of 
the null hypothesis, if in fact the RS are Poisson 
distributed. 


Fer example consider the non-homogeneous Poisson 
precess over the interval (0,100]. If this interval were 
partitioned into 100 unit segments the number of events 
occurring in each segment should be Poisson distributed. ІЁ 
51 event streams are simulated over the interval, a value 
for the test statistic d may be computed for each segment, 
noe. 


ES (n EN 
nm CE. И җе” 


n 
р 


Assuming а significance level a= .05 and comparing each 4 
Bothe critical value HE 67.5008, we would count the 
nunber cf test statistics such that q > 6/5048. I£ the 
number of events in the segments are indeed Poisson random 
variables we would expect, on the average, that 5 out cf the 


100 dp values would exceed the critical value. 


In the above discussion of the dispersion test, no 
mention was made of the parameters of the Poisson random 
Variables (assuming that they are Poisson). This is because 
the dispersion test for Poisson distributions is parameter 
independent. Thus it is necessary to perform a second 
hypothesis test to determine if the mean number of events 
per segment is equal to the difference of the integrated 
intensity function evaluated at the right and left end 
points of the segment respectively. For large k the sample 


nean n is normally distributed. Under the null hypothesis, 





a. 
i.c. that the mean of Ny EL = Up (where а, and a, are 
the end pcints of segment p) ^n is normally distributed with 


mean up and variance Up Ke The test statistic is 


P АЕА 


and the critical value is the M percentile of the 
ncrmal distribution, 44/2. Again, there would be one 
hypothesis test for each segment and if we made 1000 such 
tests we would expect to reject the rull hypothesis 


approximately 1000q times, on the average. 


Continuing the example above, it 1s now of interest 
to see if the mean number of events on each segment  agress 
With the theoretical value determined by the integrated 
intensity function. This time assuming a significance level 
По = .10, the critical value for each test statistic Co 
would be 204/27 1-655. Fach So > 1.645 would be counted and 
if in fact tte random variables have the hypothesized means, 
we would expect (on the average) that approximately 10 test 


СЕНАТЕ 5Е1Сс5 cut of 100 would exceed 1.645. 


After the dispersion test and hypothesis test for 
the means have been performed on event streams from each 
gegorithm, the test results may be compared to see if either 
algorithn appears to have null hypothesis rejection 
percentages which are not consistent with the о levels of 


the hypothesis tests. 


(ТЕ should be noted that during the development of 
the programs histograms of the event streams were plotted to 
determine whether or not they "fit" their respective 
intensity аре езй ерус Both programs seemed to produce 


Satisfactcry results by this subjective measure.) 
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СО OTHER CONSIDERATIONS 


1. Intensity Function Cateqory 


Algorithm B classifies the input intensity function 
into one of six categories and then determines what action 
to take to generate the appropriate event stream. It might 
be expected that the speed of Algorithm B will be affected 
not only by the total number of events to be generated but 
also by the Ооо ои по итеп” function. 
Decomposition of the intensity function into four sections 
will prckably require more time than if it must be divided 


into just two sections. 


АШ ОТГ1САШ treats all intensity functions alike. 
ENncretore, it is important to compare the two algorithms for 
all six possible categories of intensity functions. It is 
possible that Algorithm B could be faster than Algorithm А 
for one category of intensity function and slower than 
Кошево A fcr a different category of intensity function. 
The six intensity functions selected for the comparison 
represent all of the categories defined by Algorithm B, and 
meer in fact, those intensity functions illustrated in 
Figures 7 through 12 of Section III. These categories are 
numbered frcm I through VI and correspond directly to those 
listed in Ref. 11. 


2. Evaluation of c* Value in Rejection Routine 


The imrortance of reducing the size of the rectangle 
mio sing the intensity  £zunction A*(t) prior to beginning 
the rejection algorithm was mentioned in Section 17-р. pr 


possible c* shculd correspond to the maximum value of A*(*) 
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on the interval (0, 01 For Cases I and II (Figures 7 and 8) 
this maximum value is easily determined since it occurs at 
one of the end points of the interval. For Case III (Figure 
9) ЕЕ maximum values of AF (6) and AA (CE) Donngoecub 3t the 
point b where the interval (0,t4] has been partitioned into 
two subintervals (0,bjJ and (btol For Cases IV, V and VI, 
(Figures 10, 11 and 12) the maximum value of A*(t) is not so 
easily determined. LEWIS and SHEDLER (Ref. 11, p. 11] state 
that it is not possible to find this maximum explicitly. 
However an upper bound may be determined from the values К 
and k where A(t) < k and А (+) > k on the interval (0,*tpg ] (or 
Bubantervals (0,b] and (b, tg ] for Case VI). Then by setting 
КОО з кл 15 certain that A*(t) € c* on the interval or 
sutintervals cf interest. Figure 14 which illustrates the 
rejection routine for the intensity functions for Case V and 
VI used in the analysis, reveals that the c* computed in 
such a manner may not be very efficient as an upper bound 
Ог the rejection algorithn. One would expect some 
reduction in speed efficiency for Algorithm В when c* 


greatly exceeds the maximum value of A*(t). 


3. Designated Tolerance Level 


= a б Í EE лық ыс ш ш звис мс: жым: о ср = че — Жы оны ш» E ӘНЕС 


Algorithm A requires selection of а tolerance level 
which determines the stopping criterion for Newton-Raphson 
iterations. The speed of the algorithm could be expected to 
be very sensitive to this tolerance limit. By setting it 
too small the comparison would be prejudiced in favor of 
ШОТ Сп В. By not setting it smail enough the fidelity of 
the resulting event stream to the true intensity function 


would be impaired. 


E was determined that a tolerance level of 
-7 
ЕРНІ. ОЈ x 10 


ем егошапаке beyond its Significant digit capability in 


occasionally demanded that the сто те 
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single precision. This resulted in the failure of the 
Newton-Raphson iterations to converge, although in theory 
they would have done so if enough precision had been 
retained. A tolerance level of ELN, 1x 1079 seemed to be 
unnecessarily demanding on Algorithn A even though the 
computer could discriminate at this level. A stopping rule 
fixing Є = (Но X 107? was the compromise that insured 
reasonable accuracy of the simulated event stream without 


unfairly Lburdening Algorithm A. 
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(Acceptance region accounts for less than 
4% of the total area in the rectangle.) | 





10 20 30 40 50 


Figure l4a - Graph of Rejection Routine for the Case Y 
Intensity Function Analyzed 





10 20 30 40 


Figure 14b - Graph of Rejection Routines for the Case VI 
Intensity Function Analyzed 
(оте 14 = CIL59gSTRATION OF Rowse) TON-RecerTANGs REGIONS 
EII PRENSAS o мр CASE VEI INTENSITY FUNCTZIO:» 
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A. GENERAL 


Computer programs implementing Algorithms A and E were 
ccmpared in terms of efficiency as explained in Secticn V. 
This section documents the results of the comparison and 
discusses general observations of the author concerning the 
twc competing programs. Also recommendations for further 


study concerning Algorithm B are offered. 


DN Ло ОВЕ OF EFFECTIVENESS RESULTS 


Table I convincingly illustrates the superiority of 
Algorithm В over Algorithm A when computational speed is 
used as a measure of effectiveness. The column in Tatle I 
labeled "Time A/Time B" is the ratio of the total time 
required ty Algorithm A to produce a specified number of 
replications of the given non-homogeneous Poisson process, 
to the time required by Algorithn 3 to produce the sane 
number Сс Wreplicaticns of the process. ror all -six 
representative intensity functions the Pcisscn-decomposition 
and gap statistic algorithm is at least twice as fast as the 
time-scale transformation. For large event streams the 


Superiority cf Algorithm B is even more pronounced. 
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INE SIEOR EVENT STREAMS 


COMPUTATION 


Table: E 
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The most obvious disparity revealed in Table I is 
the difference in speed efficiency (of Algorithm E with 
respect to Algorithm A) between the group of Cases I, II and 
mers and the group of Cases IV, Y and VI. Іп the former 
group Algcrithm B was 5 to 7.6 times as fast as Algorithm A. 
For the latter group Algorithm B was only 2.2 to 2.5 times 


as fast as Algorithm A. 


This difference between the two grcups immediately 
made suspect the value c* used in the rejection routine of 
Bigorithm B. For the first group c* is set at the least 
Moper bound for \* (t). For the second group c* is not, in 
general, a least upper bound for A*(t) but is only an upper 
bcund. To determine whether or not this difference in the 
"quality" of the c* values could have such an adverse effect 
on time efficiency for Cases IV, V and VI, these cases were 
Simulated a second time with new с“ values which were 
empirically evaluated by graphical methods. These modified 
c* valuss were set close to the least upper bounds of the 
respective 1*(t) component of the intensity functions for 
each case. Marked time efficiency gains were observed (see 
Mescussicn cclumn, Table І) indicating that significant 
efficiency losses can be sustained due to the method of 
setting c* values in Cases IV, V and VI. 

‘ 

Three factors appear to influence the degree of 
speed efficiency gains of Algorithm B over Algorithm A. 
First is the selecticn of the c* value as discussed above. 
The second factor is the percentage of total events which 
EMgorrthn В must generate from the rejection routine. For 
the specific intensity functions analyzed, this percentage 
ranged from 4% in Case V to 24% in Case VI. The fewer 
events required from the rejection algorithn relative to the 
number generated from the gap statistics algorithm, the 


greater the efficiency of Algorithm B. 
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Finally the total number of events in an event 
stream affects relative efficiencies between Algorithms A 
and B. This is because as the number of events to be sorted 
ру Algorithm A grows, the less efficient its sorting routine 
becomes. Sorting requirements for Algorithm B are greatly 
reduced due to the properties of the gap statistic procedure 
used. Therefore as total events increase, so does the 
efficiency advantage of Algorithm B over Algorithn A 


increase. 


2. Computer Memory Requirements 


comm ШШЕ >; A eee = чш = A A сыш o en Vu 


Algorithm A requires approximately 72,000 bytes of 
core storage. Algorithm B demands about 92,000 bytes. Both 
programs are costly in terms of storage due to the number of 
library routines used and the need to store up to 5000 event 
times. The difference of 20,000 bytes would not ре 
important unless computer capacity were being challenged. 
Memory requirements for both algorithms could be reduced by 
reducing the event stream capacity of the programs. A 
reduction of capacity from a maximum of 5000 events to a 
maximum cf 2000 events in each program would result in 
storage requirements of 54,000 bytes and 68,000 bytes for 
Algorithm A and Algorithm B respectively. 


Be Fidelity 


MA AQ aa ee GG 


Table II lists the results of six series ОБ 
hypothesis tests of the types discussed in Section V-B. 
Three series of dispersion tests and three series of 
hypothesis tests for the mean were conducted. Intensity 
functions for Cases III, IV and VI were selected fcr these 


tests. 
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Fer each of the three cases, a large number (500, 
900 or 1000) of each type of hypothesis test was performed 
for a fixed level of significance q. The number of times 
each null hypothesis was rejected was recorded and the 
statistic à = ((nunmber of  rejections)/(total tests)} was 
computed for the series of dispersion tests and the series 
of hypothesis tests for the mean. If the null hypothesis is 


true, then a is an estimate of the significance level «а. 


From Table II it appears that for the hypotheses 
tests associated with Case III and Case IV, à gives a very 
good estimate for a, the true significance level. For Case 
VI some disparities are observed. For Algorithm B the à 
Ме for the dispersion test (i.e. hypothesis test that 
variates are Poisson distributed) is 0.039 whereas the 
significance level is а - 0.01. Algorithm A yields a much 
better @ value of 0.012. For the hypothesis tests for the 
mean, both Algorithm A and B yield a values almost twice 
maat of a. 


These results, though interesting, are inconclusive. 
For the dispersion test, the statistic 
k = 12 
E | 
су = a DA S 


n 
р 


is only approximateiy distributed as a chi-squared random 
variable. It has been shown by FISHER (Ref. 4] that for 
Poisson randcm variables with low expectations, hypothesis 
tests based on the chi-squared distribution may produce 
EDUISOUS results, Also, ке would expect that if the 
distribution of d is only approximated by the chi-squared 
distribution, a the relative error between the 
Bepeoximate and true distribution is greatest in the tails, 
i.e. at high percentile values. Case VI was tested at an a 


level c£ 0.01. It is also the intensity function with the 


Ty 





fewest expected number of events of the three intensity 
Auctions tested. Басса СЕС све interval into 
segments produced some segments which had a low expected 
number of events (see the right tail of the density function 
pictured in Figure 12). The results for Case VI rather than 
indicating any serious flaws in either of the algorithas 
suggest further research into finding better ways of 
hypothesis testing for sparse event streams and highly 


discriminating levels of significance. 
Results or the hypotheses tests on Cases III and IV 


indicate that both algorithms simulate event streams that 


ccnform to the desired hon-homogeneous Poisscn processes. 
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RESULTS OF HYPOTHESES TESTS FOR FIDELITY 


Table I1 
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С. GENERAL CBSERVATIONS 


1.  Frogramming Ease 


The programming of Algorithm A is Simpler than the 
programming of Algroithn В. Because the time-scale 
transformaticn technique handles all intensity functions 
alike, several different cases need not be identified. 
Meecorithm B must categorize the input intensity function 
into one of six categories and then must treat each category 
in a unique way. A rough indication of this difference in 
the degree cf complexity is the number of instruction 
statements required by each progran. The progran 
implementing Algorithm A required 142 instructions. The 


program fcr Algorithm B required 364 instructions. 





Algorithm B employs ап exact method in generating 
the non-hcmogeneous Poisson event stream. It is exact in 
the sense that all event times are calculated directly and 
the precision of those calculaticns are limited only by the 


Payeotcal constraints of the computer. 


Algorithm A employs a Newton-Raphson iterative 
method tc approximate event times on the interval. The 
Stcpping criterion for the technique is limited by machine 
precision considerations. Also the stopping rule does not 
give a firm account of the magnitude of error that can be 
expected in the event times. eken epsilon criterion is 


applied tc the value of the integrated intensity function 
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and not to the abscissa, (or time interval axis). For the 
integrated intensity function, F(t), 1t is conjectured that 
леп а function value F(ti;) = uj, 1f a t = t* can be found 
euch that [F(t') - (С)! SA, 13 very close to ti. 
Although this is a reasonable assumption given the 
well-behaved nature of the integrated intensity function, 
the problem cf not knowing how close t' is to t; may tend to 


reduce one's confidence in the algorithm. 


3. Initialization 


m. 


The initialization time required for the programs 
implementing the two algorithms has not yet been mentioned. 
Computation time comparisons did not include compilaticn or 
linkage times required for the two algorithms. These 
initializaticn times were significantly different, being 
approximately 16 seconds for Algorithm B and only 8 seconds 
On Algorithm A. Should simulation of only one or two event 
streams te required it is likely that Algorithm A may 
actually require less total time than Algorithm B. The 
difference between the initialization times could probably 
te reduced if the programs were to be rewritten in Assembler 


language. 


D. CONCLUSICN AND RECOMMENDATIONS 


ШЕН Conclusion 


For general usage, the Poisson-decomposition and gap 
Beatistic method of LEWIS and SHEDLER [Ref. 11] is clearly 
superior to the time-scale transformation methcd in 


generating a non-homogeneous Poisson process with a 
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degree-twc exponential polynomial intensity function. 


The main advantages of the Poisson-decomposition and 
gap statistic algorithm are its speed and its qualification 


as an exact method of variate generation. 


The time-scale transformation method enjoys some 
advantage in core storage requirements and in lower program 
initialization time. These advantages are not sufficient to 
overcome the relative deficiencies of much greater 
computation time and uncertainty about the accuracy of the 


approximating mechanism for determining event times. 


D. Recommendations 


The full potential of the Poisson-decomposition and 
EU statistic algorithm can not be realized until it 
includes a better method for selecting an upper limit value 
C* for Cases IV, Y and VI intensity functions. Algorithms 
for choosing a c* value that wiil be close to the maximum 
value for the function A*(t) should be developed and 
Encorporated into Algorithm B. A single variable search 
technique (such as a Golden Section search or Bisection 
search) for finding the maximum value of a unimodal function 


may be appropriate. 


The ccmputer program for Algorithn B should be 
rewritten in Assenbler language in order tc reduce progran 
initialization tine. The progran could then be developed 


into a library routine for general use. 


The question ОТ fidelity of the simulated 
ncn-homogenecus Poisson process to the true theoretical 
process should be investigated further. Of special interest 


is the question of how well the algorithm simulates sparse 
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event streams. Methods for conducting the dispersion test 
and the hypothesis test for the mean at very high levels of 
Significance (i.e. a = .01, .005) for intervals with a low 


mean number of events are needed. 


83 





APPENDIX A 


БОР ОЕ ТПЕ NALEDITY ОР SCALING THE INVERSE DISTRIEUTION 
FUNCTION 


Propositicn: Let T be a continuous random variable with 


UP FEO <P EN cm oo 


ÍmMetrabution function Ao Such стас 
T 


Flt) = 0 

ET 000 
F(t) = IA mec ER 
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Let Y = ror pf?) such that T кол ) and 15 а real number. 
Then the random variable Y is unifornly distributed cn the 
interval (0,50) - 
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0) then 
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Therefore, 


ES D SMS 


and Y is uniformly distributed on the interval (0, ug) - 


85 





APPENDIX B 


ERROR FUNCTICN AND DAWSON'S INTEGRAL REPRESENTATION CF THE 
INTEGRATED INTENSITY FUNCTION A(t) 


te Standard form for the error function is: 


€ -u* 
= du 


sl" 
O 


Dawson's integral is usually represented as: 


-2 E u? 
ДЕ е qu 


Both of these integrals may be calculated to any desired 
degree СЕ accuracy by using the exponential series 


expansion, 


X 
lag 
s 


u u | 
Thus, е апа е may be written as 
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and 


ú Ш 
e du = » ram c 
= b n 
n= 
со ¿2n+l 
А E (2n+1)n! 


8 = 2 , 
Mubtiolyicg by t results in 


С 2 co 2п-1 
=2 u t 
0 O (Ana 


and the series representation for Dawson's integral is 
revealed. This series will converge for all t, althouch the 


value of the integral increases very rapidly as t increases. 
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Using the same argument it is easily shown that the series 


Р t  -ue. : 
representaticn for Lo e "au is 


со 


(1) 0 „en+l 
(ара! — ^ 
EO (2n+1)n! 


The error function 15 obtained by multiplying the integral 
ЕЕ Бе ccnstant 2//m  ,i.e., 


t NC со 
| е ШІП - | 
0 = 


|x 


Although the series expansion may be used to calculate both 
functions, there are very efficient computer library 
routines available for computing Dawson's integral and the 


БЕБОК functicn. 


The FORTRAN supplied procedures ERF and DERF [Ref. 15] 


evaluate the above function for input values of t. 


Ше ІМ5І (International Mathematical and Statistical 
Libraries, Inc.) FUNCTION MMDAW (Ref. 6] evaluates Dawson's 


integral for input values of t. 


ESE remains to be shown that the intensity function 


i 2 
à (t) = exp (0, + at + a, 


interval (0,t,] using either ERF (or DERF) or MMDAW. 


t} may be integrated over the 


88 





Consider, 


0 0 


2 
1105540) — | ME ! exp(ay + ayt + at )at 


By completing the square of the exponent the expression 


beccmes 


2 
Eo с = 
A (tg) = || SE ae | ехр | aA (t + p Gig 


(A(0) 0, since 0 is the left end point of the interval 


over which 1(t) is defined). 








For a. > 0, 
a % | ie 
Alt,) = exp ya, - ia, f К C dt 
2 
МЕ ing, 
лыш де + =, 
2 20.5 
du = va, dt 
a 
=> - n 
K = exp E e 
and adjusting the limits of integration, the expression 
becomes: 
b 2 b 2 a 2 
A(tg) = 2 f e” du = = | f e au- f е“ a 
vo, a Ya, 0 O 





where, 


and 
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A(t.) = E (b*b ^? f e` du- deae f e” du} 


0 Ie, 0 0 


Inside the brackets are two Dawson integral expressions 
multiplied by constants. Outside the brackets is just one 
more easily determined constant. Using  MMDAW twice and 
multiplying by the appropriate constants will give the 
desired integral value. 


ШЕ о. < 0, 


2 
af] “o a, 2 
A (tg) - exp bn | ехр = е5 T ee: e 
à 0 2 
at EQ a, 2 
е к | y E 20,1 dt 
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Substituting as before gives 
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and the errcr function can be recognized within the 


brackets. 


Two error function calculations, a subtraction and a 


multiplication by a constant are sufficient to determine 


е). 
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APPENDIX C 


COMFARISCN OF GAP STATISTICS ALGORITHM AND TIME-SCALE 
TRANSFORMATION ALGORITHM FOR SIMULATION OF NON-HOMOGENEOUS 
POESSONSEROCHRSSES WITH LOG-LINEAR INTENSITY FUNCTIONS 


LEWIS anā  SHEDLER [Ref. 13] present two algorithns 
for simulation of non-homogeneous Poisson processes with 
intensity function \(t) = exp(ß, *B,t)- The first algorithm, 
which uses a time-scale transformation is easy to employ 
since the inverse of the integrated intensity function is 
known explicity. The second algorithm uses gap statistics 


to generate event streams. 


Both methods are exact. That is, cther than the 
physical precision constraints inherent to the computer, no 
approximations are necessary in generating event times. 
However the gap statistics algorithm obviates the need for 
ordering an array of variates. Also, the gap statistics 
algorithm takes advantage of a fast standard exponential 
Variate generator, (Random Number Generator Package 
LLRANDOM, Ref. 12), thereby avoiding the time consuming 
ccmputation cf taking logarithms. The time-scale algorithm 
must perform both the ordering of variates and the 
computation cf logarithms. 


Bcth algorithms were programmed and tested for 
computaticnal speed efficiency. The gap statistics 
algorithm was roughly twice as fast as the time-scale 


eeenstormaticn algorithn. 
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The SUBROUTINE NHPP2 uses the gap statistics 
algorithm within the Poisson-decomposition and gap statistic 
algorithm (Algorithm B). It is the presence of SUBRCUTINE 
NHPP2 which markedly reduces the ordering requirements of 


Algorithm B from those necessary in Algorithm A. 
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APPENDIX D 


POISSON VARIATE GENERATION 


A. | BACKGROURÓL 


Acknowledgement: The content of this appendix is a 


paraphrase cf portions of Chapter 11 of an unpublished book 
by J. Н. AHRENS and 0. DIETER [Ref. 1]. 





The Pcisson distribution has the probability mass 
Function 
= 
= e А 
Pi i! 
where р; is the probability that exactly 1 events are 


observed to cccur ina unit time interval, given that events 


arrive at an average rate A. 


The intervals between events in a Poisson process of 
rate д) are independent and have an exponential distrikution 


meen Mean 1/% i.e. 


f(x) = ле 


The prokatility integral transform method can be used to 
obtain a sample of exponential variates from a sample of 


uniform (0,1) variates, Ur Ur > ‚ »« , Since 





A A rate event stream can then be simulated from the 


exponential variates using the following interpretation: 


15+ event occurs at En 
2nd event occurs at i: * x 


etc. 


This property yields two simple methods for generating a 
realizaticn cn the Poisson random variable. These are the 


additive method and the multiplicative method. 


ПРЕ ADDITIVE METHOD FOR GENERATING A POISSON VARIATE 


The number k of unitS arriving in a unit interval ina 
Poisson process with rate A must be found. This will be 


the k for which 


ЕТ СТ 


JE 2 К 
апа 
x] Bo quur иа i 
or equivalently, for which 
-2n а, - 2n U, -...- Qn чу A 
and 
-2n uy 7 Qn СЕ = gn үз Qn Uy 11 > x (1) 


(Ncte: Since u, < Iyi 222035 4. > 0) 


Thus by using uniform variates, computing logarithms, 


adding and ccunting, the Poisson variate is obtained. The 
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problem with the additive method is that if \ is very 
large, say 2000, it requires considerable effort to generate 
a single Poisson variate. Although addition is extremely 
fast on a computer, the computing of logarithms is 


relatively slow. 


ШЕРТПЕ MULTEPLICATIV: METHOD FOR GENERATING A POISSON 
VARIATE 


ми етапа (1) by -1 gives; 


Un lu, u, © ео о e uy) 2 -À 
and 
An (u, "u, o • • • • Uy tU. 1) < -À ; 
or equivalently 
u • 11 • e ө o o u > Sas 
1 е = 

апа = 

uz "u, • ео ө ee = 5 


501 by generating uniforms, successive multiplication and 
eounting, the Poisson variate NS | protec ede The 
multiplicative method eliminates tae need for computing 
logarithas. However +гоцріз is encountered for large 
femmes since © becomes very small аз А increases. 
ШОССЕ сох "problems occur :or A values of approximately 175 
ата Тагдег Оп the IBM 360/67 computer system using Single 
precision word lengths. The precision problem may De 


Overcome by partitioning the Poisson random variable into 


the sum of two or more Poisson random variables. However 
@memaverage number of multiplications and uniform variates 
шесиштед 1S Still proportional to A. More erficient 
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methods exist such as the Gamma method. 


D. THE GAMMA METHOD FOR GENERATING A POISSON VARIATE 


(The fcllowing discussion is a condensed account of Ahrens 


and Dieter's presentation; Ref. 1, pp.11-20,21, 22) 


In order to obtain a sample k from the Pcisson 
distribution of mean U „ select a positive integer n 
Meyetcaliy n is a little smaller than U). Then take a 


sample x from the Gamma distribution of parameter n. 


len X OU return a 
sample k from the binomial 
distribution with parameters 
п-1, Н /х 


(С) TE х < р take a sample 
J from the Poisson 
distribution of mean u-x and 


cecuri K= nt j. 


The sample x simulates the n-th event (arrival) ina 
Poisson piecess of rate 1. If x >u then there are n- 
arrivals in the interval (0,x), and each of these has a 
probability cf u/x of being below Ы [Case (1) 1. TE 

ms then the n simulated arrivals are all before U, 
and the sample j indicates the additional events between х 
and uy [Case (2) ]. 


(See thie point Ref. 1 includes a formal proof of the 


procedure, which will ke omitted.) 
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The explicit algorithm relies on an efficient methcd for 
sampling frcm the Gamma distribution. (Such a method has 
been implemented as W. R. CHURCH Computer Center  litrary 
routine GAMA by D. W. Robinson and P. А.Н. Lewis and is 
documented in Reference 10.) The constant n is 
arbitrarily сһоѕеп аѕ п = [0.875 р]. This avoids the Case 
(1) almost ccmpletely if u is large and a simple method for 
sampling from the binomial distribution becomes 
sufficient: the Bernoulli process is simulated in Steps 
8-10 belcw. The algorithm for sampling from the Poisson 
distribution with mean u- x in Case (2) is the simple 
multiplicative method explained above andis employed in 
Steps 3-5. However, if u - x is still large (larger than 
Lnes'"cut-cff point" c) the entire procedure is iterated with 
a new auxiliary sample x from a Gamma distribution of mean 
ПИ (0.875 (U ~- х) 1] etc. 


Algorithm - The Gamma Method, Ahrens and Diete 


1. Initialize k + 0, W + у. 
2 PA OO 6 (c = 15 Was’ determined 
empirically Ey Ahrens and Dieter to be 


reasonable). 


з (Start Case (2)). Set op + 1 and calculate 


-W 
Б <= є А 


4. Generate u and set p * peu. If p<b deliver 


98 





I: Increase k «+ K+1 and go to 4. 


б Set n -[dw] where d - 7/8 (d value was also 
selected empirically). Take a sample x from the 


Gamma (n,1) distribution. If x» w go to 8. 


p Set k = k * n, v — w-x and go to 2. 


9% (Start Case (1)). Set p — w/x. 


9. Generate u. If u$ p increase k -k * 1. 


Мени n=. ІР п > 1 go to 9. 


11. Deliver k. 


Ahrens and Dieter claim that the computation time for 
the Gamma method grows ultimately like а + в 1а у. 
Computaticn time for the additive and multiplicative methods 
grcw like U . Therefore for the large u values typically 
demanded in the simulation of non-homogeneous Pcisson 
processes, the Gamma method is signally superior in terms of 


ccmputaticn speed. 
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APPENDIX E 


GRAEHICAL PRESENTATICN OF NEWTON-RAPHSON METHOD 


A. GENERAL PRINCIPLES 


~ = = amm = dub ~ = ~ 4. 


t; (unknown) t 


Figure E-] 


швадјег the function F(t) pictured above in Figure E-1. 


Memon jective is to find, for a known value wu; оп the 
vertical axis, a unique corresponding value t; such that 
F(ti) » uj. If an explicit expression for the inverse of 


F(e) exists, the calculation may be made directly, i.e 


ө „ә э Ff 
t, = Fl). Otherwise, numerical methods must be used to 
L 1 
EL roximate t. 
l 
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Figure E-2 
Form the new function G; (t ) = F(t) - u, (see Figure 
Е-2). In effect, the horizontal axis has been disrlaced 


upward a distance u;. Now if a t* can be found such that 
G; (t*) = 0, then F(t*) = 0. and t* = tir the desired value. 


Assume that an initial approximation for t* can be made, 


say t! 4 If Gi (t ) and Gi (ti ) = g; (ti ) can be computed 


we can write 


G iti) 
q GER) = г r , 
t (51 t5) 
where tj is the intercept of the line tangent to aree) at 


G; (t? ) with the t-axis. Solving the above expression for 


gives: 
2 


It is evident from the graph that © is closer than е 
КО Те unknown t*. 








Figure E-3 


Repeating the procedure (Figure E-3), using tj as a new 
merroximaticn will yield t3 a still Detter approximation of 
a, in general, given an approximation ty has been 
obtained, a closer approximation ty} сап be calculated by 


the relationship 


It can be shown that successive approximations will converge 
to the value t* provided G; (t) satisfies certain criteria 
mist. 5, р. 449,450]. 


Since Gi (t) Hohn onm ol a distribution function F(t) 
it satisfies ali the conditions for convergence except for 
the case where the interval [ti ac contains an inflection 
punt. If the inflection point(s) can be identified, this 
troublesome situation can easily be avoided, by handling the 
function G, (°) in discrete sections over the discrete 


intervals, (0,t* ЖА (t? ‚ty Mog ET (V yvtk ) where the El 
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Пи паге values suck that G;(t* ) are all the 
inflecticn pcints of Gi(e) over (0,50). (For the special 
case under consideration at most one inflection point will 
be encountered. It is the maximun or minimum of the 


intensity function.) 


Successive approximations are computed until 
БЕРДІ << at which time E is assumed to be close 


encugh to t*. 


B. INITIAL APPROXIMATIONS 


The prceklem of obtaining a good initial approximation 
for each t; was solved in the foliowing manner (see Figure 
E-4 below): 


= = -- = = => <= =з <> <= => ED A <> = = <> => == => = aa «a == <> woh 


= 
Co 
= 
fa 
4 
сл 
ct 


22 


Figure E-4 
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The interval (0,tg] was divided into several segments 
(Ort, Je (т ет], etc. Then the values 5 (3 = 1,27...) 
were computed. Each resulting abscissa and ordinate value 


was stored in an array, i.e., 





For each value u the array «sas searched until two function 


values were located such that Е (2) TS РТ) + Either 


J 
E ог ETSI was then selected as the initial approximation 
for the value ti (i.e. ty 1) such that Ee.) - u;. The 
function G; (t = кав then  £crmed and the 
Newton-Raghscn iterations performed until the epsilon 


Ci iterion was satisfied. 
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